
• r 

; r : : 





Manuel D. Salas 


//l/ ~D2~~ 

2 , 035 ^ 


f. } Vi 


Shock Wave 
Interaction 
With an Abrupt 
Area Change 


(NASA-TP-^113 j SHOCK WAV? INTERACTION WITH 
AN ABRUPT AREA CHANGE (NASA) 16 p CSCL 01A 



N9 l - 1 7 i.40 



4 


I 


■ _ f . h y - ' K. .. , 


.1 • ' — : 


1 ' — 


Unci as 

Hl/02 0030859 







NASA 

Technical 

Paper 

3113 


1991 


Shock Wave 
Interaction 
With an Abrupt 
Area Change 


Manuel D. Salas 
Langley Research Center 
Hampton, Virginia 


NASA 

National Aeronautics and 
Space Administration 

Office of Management 

Scientific and Technical 
Information Program 



Abstract 


The wave patterns that occur when a shock wave interacts with an 
abrupt area change are analyzed in terms of the incident shock wave 
Mach number and area-jump ratio. The solutions predicted by a self- 
similar model are in good agreement with those obtained numerically 
from the quasi- one- dimensional time- dependent Euler equations. The 
entropy production for the wave system is defined and the principle of 
minimum entropy production is used to resolve a nonuniqueness problem 
of the self-similar model 


Introduction 

The interaction of a shock wave with a channel of 
rapidly varying cross-sectional area is of interest in a 
number of practical problems, such as the passage of 
shocks through wire-mesh screens, the starting pro- 
cess in a supersonic wind tunnel, and the phenomena 
that occur in piston engines and jet engines. Previous 
investigators (refs. 1-3) have shown that a self-similar 
inviscid model with a discontinuous area change can 
provide good agreement with experimental observa- 
tions. A solution to this model is obtained by guess- 
ing a self-similar wave pattern with its origin at the 
location of the area discontinuity. The guessed pat- 
tern is validated if the conservation laws of mass, mo- 
mentum, and energy can be satisfied. The problem 
essentially depends on two parameters: the strength 
of the incident shock wave, measured by the shock 
wave Mach number, and the area ratio across the 
discontinuity. This parameter space is rich in the 
number of possible wave patterns and several inves- 
tigators (refs. 3-5) have indicated that more than one 
wave pattern might satisfy all the conservation laws. 

The existence of multiple solutions was the sub- 
ject of an article by Oppenheim, Urtiew, and Stern 
(ref. 5). They showed that, in a region of the pa- 
rameter space corresponding to supersonic flow be- 
hind the incident shock and within a certain range 
of area contraction, three wave patterns could sat- 
isfy all the conservation laws. Oppenheim, Urtiew, 
and Stern conjectured that the ambiguity could be 
resolved by invoking the minimum entropy produc- 
tion principle. This led them to accept two solutions 
in this region, one with a standing shock wave within 
the area contraction. Rudinger (refs. 6 and 7) ques- 
tioned their conclusion, pointing to the well known 
fact that a standing shock in a converging channel 
is unstable. Through a study of the transient phe- 
nomena produced by a steep, but continuous, area 
variation, Rudinger concluded that the only solution 
that could be realized in this ambiguous region cor- 
responds to a wave pattern with a rarefaction swept 
downstream. Here we show that the solution pro- 


posed by Rudinger can be reconciled with the mini- 
mum entropy production principle if the entropy pro- 
duction is properly defined. 

Rudinger 5 s transient analysis was based on a 
graphical method of characteristics. This tedious ap- 
proach limited Rudinger to the study of three specific 
examples. In order to establish conclusively that a 
reflected shock wave cannot be formed in the region 
of ambiguity, Rudinger proceeded to show that the 
waves reflected from the transmitted shock cannot 
coalesce until the head of the reflected wave becomes 
stationary, that is, the flow becomes sonic. Implicit 
in the proof is the assumption that the head of the re- 
flected wave becomes stationary for conditions on one 
of the boundaries of the region of ambiguity. While 
this is true for the self-similar model, it is not clear 
that this is also true for the transient problem. 

For an area divergence, no multiple solutions are 
known to exist. The region of ambiguity that occurs 
for an area contraction can be shown to extend into 
the region corresponding to an area divergence in 
parameter space. Here, however, a unique solution 
with a standing shock is found. 

The purpose of this paper is to map the differ- 
ent wave patterns that take place for the self-similar 
model in terms of the incident shock strength and 
area ratio and to verify the validity of these solutions 
by solving the time-dependent quasi-one-dimensional 
Euler equations for flow in a channel with a steep 
cross section. The study is limited to monotonically 
increasing or decreasing areas. The problem is de- 
fined and its method of solution is explained in the 
first section. This section also investigates the flow 
patterns that take place for an area divergence and 
an area contraction. Following in the next section, 
the quasi-one-dimensional model is introduced and 
the numerical method for solving this problem is out- 
lined. The results section compares the self-similar 
model and the quasi-one-dimensional model. Finally, 
conclusions are discussed in the last section. 


Symbols 


A 

channel area 

Al 

channel area to the left of 
area discontinuity 

Ar 

channel area to the right of 
area discontinuity 

a 

speed of sound 

C 

Riemann variable defined by 
equation (21) 

D 

Jacobian matrix defined by 
equation (27) 

e 

specific total energy 

F 

flux matrix defined by 
equation (24) 

K 

constant in minmod limiter 

L(/ 

left eigenvector matrix of C 

M 

Mach number 

M* 

value of Mi corresponding 
to sonic conditions in 
region 3 

P 

pressure 

Pr 

pressure ratio (see eq. (9)) 

Q 

source vector defined by 
equation (24) 

n 

residual defined by equa- 
tion (33) 

s 

entropy (see eq. (5)) 

t 

time 

u 

unknown vector defined by 
equation (24) 

u 

velocity 

W 

characteristic variable 
vector defined by equa- 
tion (28) 

w 

value of W returned by 
minmod limiter 

w 

shock wave speed 

X 

axial coordinate 

y 

argument for minmod 
limiter 


z 

argument for minmod 
limiter 

a 

area ratio (see eq. (1)) 

ot c 

asymptote of curve c (see 
%• 2) 

<*d 

asymptote of curve d (see 

fig- 2) 

0 

defined by equation (25) 

1 

specific heat ratio 

6 

defined by equation (4) 

K 

defined by equation (4) 

A 

characteristic slope 

V 

constant appearing in 
minmod limiter 

p 

density 

a 

constant appearing in 
equation (22) 

X 

entropy production 

Subscripts: 

i 

incident shock 

k 

time counter 

n 

space counter 

r 

reflected shock 

t 

transmitted shock; differen- 
tiation with respect to time 

X 

differentiation with respect 
to X 

0 

starting conditions 

1,2,3,... 

regions of flow 

Superscripts: 

k 

time counter 

+ 

forward difference 

- 

backward difference 

0 

Runge-Kutta stage 

Special notation: 

a, b, . . e 

curves in figure 2 

I, II, . . IV 

quadrants in figure 2 

la, lb, Ic, Ha, lib, 
Ilia, 111b, IVa 

flow patterns in quadrants 
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Self-Similar Model 

Consider two infinitely long constant area ducts 
that are connected by a short, monotonically increas- 
ing or decreasing transition section. Assume that the 
transition section is small enough that it can be re- 
placed by an abrupt transition. Further assume that 
the gas inside the duct is at rest. We are interested 
in establishing the valid wave patterns that result 
when a shock wave moving from left to right passes 
through the discontinuous area change. Let x = 0 
be the location of the area jump, and let t = 0 be 
the time at which the incident shock reaches the area 
jump. Because there is no reference length, we expect 
the solution to be constant along rays originating at 
(0,0). That is, the dependent variables are only func- 
tions of the ratio x/t. 

Method of Solution 



i 

1 

I “ 

Figure 1 . Typical wave diagram for the interaction of a shock 
with an area discontinuity. 


A typical wave diagram of the interaction of a 
shock wave with an area discontinuity is shown in 
figure 1. In all such figures that follow, the area dis- 
continuity is depicted as a long-dash line, the shock 
waves are depicted as thick solid lines, a contact sur- 
face is depicted as a short- dash line, and an expansion 
fan is depicted by thin solid lines. Region 1 is the re- 
gion to the right of the area discontinuity and ahead 
of the transmitted shock; region 2 is the region to the 
left of the area discontinuity and ahead of the inci- 
dent shock. The flow is assumed to be at rest in both 
of these regions, and the pressure and density are as- 
sumed to be uniform. The pattern shown in figure 1 
is one of many that we will be discussing later. The 
flow conditions leading to this pattern correspond to 
a high incident shock Mach number and a high area 
ratio. The area ratio a is defined as 



(i) 


where Ai is the area to the left and Ar is the area 
to the right, both assumed to have a nondimensional 
length of 1. 

The conditions in region 3, immediately behind 
the incident shock, are evaluated from the Rankine- 
Hugoniot relations: 


Here, u, p, a, and p are the velocity, density, speed 
of sound, and pressure, respectively. Pressure and 
density are nondimensionalized by their initial values 
in region 1, and all velocities are nondimensionalized 
by the speed of sound in region 1 divided by y/j. The 
subscripts in equations (2) denote the appropriate 
region. The Mach number of the incident shock is 
denoted by M h and is given by 



(3) 


where tCj is the incident shock speed. In the follow- 
ing, 6 and n are given by 


6 = 


7~ 1 
2 

7+1 


(4) 


From the definitions of the speed of sound and the 
entropy, we have 


Fm 

a 3 = \ 

V PS > 

S 3 = ln(p 3 ) - 7ln(p 3 ) . 
Conditions in regions 1 and 2 are given by 


(5) 


a 2 


n 3 = 


K - Q 1 


K.Mi 

P2KM? 


p 3 05A/2 + 1) 

yM ? - 6 

P3 = P2 


(2) 


Ul 

= u 2 

= 0 

ai 

= a 2 

= V7 

P2 

= Pl 

- 1 

P2 

= Pi 

= 1 

s 2 

= Si 

= 0 
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The flow in region 3 becomes sonic when M l 
equals some critical value M* . If we set M3 = 1, 
using equations (2) and (5), we get 


, f *2_ ( 7 -7) + V( 7 -7) 2 - 16 ( 2 -7) 

* 4(2- 7 ) 


(7) 


Prom the conservation of total enthalpy, 


<*6 — °7 



1/2 


and since the flow is sonic in region 6, 


( 12 ) 


For values of M x greater than A/* , the flow in region 3 
is supersonic. For 7 = 1.4, M? = 2.068. As Af* — > 00, 
the Mach number in region 3 approaches the value 
I/V7?. For 7 = 1.4, the upper limit for M3 is 1.890. 

Across the contact surface, the following two re- 
lations must be satisfied: 

P4 = P5 I 

(8) 

U4 = «5 J 


If the Mach number Mt of the transmitted shock 
is known, then the flow in region 4 is defined by 
equations (2) and (5), with A/j replaced by Mt and 
subscripts 2 and 3 replaced by 1 and 4, respectively, 
The Mach number of the transmitted shock, in terms 
of the pressure ratio p r = p\jp \ , is given by 


uq = a 6 (13) 

The density and pressure follow from equations (5): 

ln { a lh) ~ s 6 


PC, = exp 


26 


P 6 


_ P&4 


(14) 


The Riemann variable on the characteristic with 
slope u + a, crossing the expansion fan, provides one 
piece of information about region 5. If we guess the 
slope of the expansion tail, A5 = — <15, after some 

simplification we get 


M t = 


I Kp r -f 5 
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(9) 


Therefore, with £>5 known, region 4 is completely 
defined. 

In general, the wave pattern between regions 3 
and 5 will be different from that shown in figure 1. 
The specific pattern will depend on the value of the 
incident Mach number and the area ratio. Here we 
illustrate how the solution for the wave pattern of 
figure 1 is obtained, with the understanding that 
similar procedures are used as the wave pattern 
changes between regions 3 and 5. 

The Mach number in region 6 (region 6 is actually 
one point in space), immediately to the right of the 
area discontinuity, is sonic. Therefore, by solving the 
conservation of mass relation written in the form 


M 6 ( i+ 6M* \ k/2B 
M 7 \ 1 + 5 a/2 J 


( 10 ) 


we can obtain A/7. Given A/7, the Rankine-Hugoniot 
relations across the reflected shock can be solved it- 
eratively to obtain the solution for region 7. With 
region 7 defined, we turn our attention again to re- 
gion 6. Since the flow is isentropic between regions 7 
and 6, we have 

S 6 = S 7 (11) 


^5 = a 6 + 


K 


a o — w 5 ^5 


(15) 


Because S5 = Sg, the pressure and density in region 5 
can be obtained from equations (14) with an appro- 
priate change of subscripts. If ^5 matches 114 the 
problem is solved. Otherwise, we continue iterating 
on A5 until U5 = 114. 

The lines — 2.068 and a = 1 lead to a natural 
breakup of the parameter space Mj, a into four 
quadrants, as shown in figure 2. In the following 
two sections, we explore the various wave patterns 
that represent solutions in each of these quadrants. 


Area Divergence 

Consider the first quadrant, Mi < 2.068 and 
a < 1. For weak incident shocks, a weak rarefaction 
wave is reflected when the shock crosses the area 
discontinuity. The effect of the rarefaction is to 

accelerate the flow before it enters into the area 

divergence. Because the flow remains subsonic as 
it reaches the area divergence, it is decelerated as 

it crosses into the big chamber. In general the 

transmitted shock is weaker than the incident shock. 
Figure 3 shows the wave pattern that is valid in this 
region, which we label la. The flow conditions for 
this figure are M x = 1.100 and a = 0.5. 
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As the strength of the incident shock increases, 
the rarefaction wave becomes stronger, eventually 
creating sonic conditions at the entrance to the area 
divergence. The locus of points corresponding to 
sonic conditions at the entrance to the divergence 
is shown as curve a in figure 2. The wave pattern 
along this curve is of type la. Figure 4 shows the 
pattern for M \ = 1.303 and a = 0.5. If a — ► 0, 
curve a approaches asymptotically a value of 1.154 
for 7 = 1.4. 

If the shock strength continues to increase, a 
standing shock develops where the area jumps. If 
we model the area change by a continuous variation, 



Figure 4. Wave pattern along curve a, type la. Mi = 1.303; 
a = 0.5. 



then as the incident shock strength increases, the 
standing shock becomes stronger and moves from the 
entrance of the divergence, where the area is A £, to 
the exit, where the area is Ar. If the area change 
is modeled by a discontinuity, the standing shock 
has no distance to move as the incident shock gains 
strength. The shock motion can only be accounted 
for through a change in the Mach number ahead of 
the shock. This in effect models the shock motion be- 
tween Al and Ar. Figure 5 shows wave pattern lb 
corresponding to a standing shock wave. The condi- 
tions for this case are Mi ~ 1.500 and a = 0.5. The 
Mach number immediately ahead of the divergence 
is sonic. From sonic conditions, the flow is isentropi- 
cally accelerated to Mach 1.927, corresponding to an 
area ratio of 0.629. After the flow crosses the stand- 
ing shock, the Mach number becomes 0.591. The flow 
is then isentropically compressed to Mach 0.427, cor- 
responding to an area ratio of 0.794. This completes 
the overall area divergence ratio of 0.5. A discontin- 
uous area change causes a squeeze of all these Mach 
number jumps into one point in space. 

As the shock strength continues to increase, the 
standing shock reaches the exit of the area diver- 
gence. At this point, the flow in front of the shock 
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.522 




goes through an isentropic expansion corresponding 
to the full area jump. The locus of points correspond- 
ing to this condition maps to curve b in figure 2. 

The wave pattern changes to type Ic with a 
further increase in shock strength. Now the standing 
shock is swept downstream, the result being the 
pattern shown in figure 6 for = 1.850 and a = 0.5. 
This pattern occurs in the region bounded by curve b 
and line Mj = M* . Above curve 5, the flow entering 
the big chamber is supersonic. As M{ approaches M* 
the reflected expansion fan disappears. 

Consider the second quadrant, Mj > 2.068 and 
a < 1. In this quadrant the flow behind the incident 
shock is supersonic. For area ratios to the right of 
curve b the pattern that occurs is shown in figure 7. 
The figure is drawn for M t = 2.500 and a = 0.5. The 
significant features of this pattern, labeled Ila, are 
the absence of a reflected wave and the appearance of 
a downstream running secondary shock. As discussed 
previously, the Mach number behind the incident 
shock is bounded by the value 1.890 for 7 = 1.4. This 
Mach number limitation does not apply to the flow to 
the right of the area divergence. Here very high Mach 
numbers can be achieved by decreasing the area ratio 
a, but keeping it to the right of curve b . For example, 
for the conditions of figure 7, Mach 2.230 is achieved 
in the big chamber. If the area ratio for this case 
is lowered to 0.15, Mach 3.512 is achieved in the big 




Figure 9. Wave pattern along curve c, type Ilia. M l — 3.500; 

a = 1.157. 

chamber. This fact was used by Ilertzberg (ref. 8) to 
design a new shock tube for hypersonic flows. If, at 
a given the area ratio is less than or equal to the 
ratio corresponding to curve 5, then the secondary 
shock becomes a standing shock. This pattern is 
labeled lib. It is very similar to pattern lb, figure 5, 
except that the flow behind the incident shock is 
supersonic and there is no reflected rarefaction wave. 

Area Contraction 

Consider the third quadrant, Mj > 2.068 and 
a > 1. If the area ratio is large, wave pattern Ilia 
occurs. This is illustrated in figure 8 for Mj — 3.500 
and a = 1.3. In this region the reflected wave is 
a shock. The subsonic flow behind the reflected 
shock is accelerated to sonic conditions by the area 
convergence. The flow is then further accelerated by 
a rarefaction wave running downstream. In general, 
the transmitted shock is stronger than the incident 
shock. If we decrease the area ratio, holding Mj 
fixed, we reach curve c of figure 2 when a — 1.157. 
The wave pattern at these conditions is illustrated 
in figure 9. It is clearly a type Ilia pattern. If we 
continue to decrease the area ratio, holding M ly we 
reach curve d when a = 1.086. At these conditions 
the reflected shock becomes a standing shock, which 
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Figure 11. Wave pattern along curve d, type Illb. = 

3.500; a = 1.086. 

is the limiting case of pattern Ilia. Curve d consists 
of the locus of points for which the reflected shock 
becomes a standing shock. Oppenheim, Urtiew, and 
Stern (ref. 5) showed that as Mj — ► oo, curve d 
approaches an area ratio given by 

a d = ^=[ 1 {\-6))- l l 26 (16) 

For 7 = 1.4, takes on the value 1.543. 

If the area ratio is just slightly greater than one, 
then we have a type Illb wave pattern. This wave 
pattern is illustrated in figure 10 for M* = 3.500 
and a = 1.06. Under these conditions, the flow 
reaches the area jump at supersonic speed. The 
area contraction compresses the flow isentropically, 
but not sufficiently to make the flow subsonic. Once 
within the small chamber, the flow is accelerated by a 
rarefaction wave running downstream. If we hold M x 
fixed and increase the area ratio, we reach curve d 
when a — 1.086. The wave pattern is illustrated 
in figure 11 and is clearly a type Illb pattern. If 
we further increase the area ratio, we reach curve c 
when a = 1.157. At these conditions, the area 
ratio isentropically compresses the supersonic flow 
behind the incident shock to sonic conditions. Thus, 
the head of the expansion running downstream in 



Figure 12. Wave pattern along curve c. type Illb. M* = 3.500; 

a = 1.157. 

the small chamber is sonic. This wave pattern is 
illustrated in figure 12. Curve c represents the locus 
of points for which the area ratio produces sonic 
conditions after the area jump. Oppenheim, Urtiew, 
and Stern (ref. 5) also showed that as M x — > oo, 
curve c approaches an area ratio a c given by 

«.= v^(|)' 1/2S (17) 

For 7 — 1.4, a c takes on the value 1.193. 

The region between curves c and d is the region 
of ambiguity discussed by Oppenheim, Urtiew, and 
Stern (ref. 5) and Rudinger (refs. 6 and 7). As we 
have already seen, wave patterns Ilia and Illb coexist 
in this region. In addition, a third pattern with a 
standing shock within the area contraction and an 
expansion running downstream is also a solution of 
the self-similar model. 

Oppenheim, Urtiew, and Stern (ref. 5) invoked 
the principle of minimum entropy production to re- 
solve the ambiguity. For each solution in this region 
they defined the entropy production \ to be 

X = ma x(S 4 ,5' 5 ) (18) 

The resulting entropy production is shown in fig- 
ure 13 for M t — 3.500 and 1.086 < a < 1.157, From 
this, they concluded that wave pattern Illb was valid 
for area ratios slightly greater than those on curve d. 
However, at some point within the region of ambigu- 
ity the standing shock pattern would take over until 
curve c was reached. Rudinger (ref. 6) objected to 
their conclusion, dismissing outright the minimum 
entropy principle and correctly pointing out that for 
an area contraction a standing shock solution is un- 
stable, as has also been shown in other investigations 
(ref. 9). Rudinger (refs. 6 and 7) further showed that 
if the area discontinuity is replaced by a steep area 
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Figure 13. Entropy production in region of ambiguity. (Based 
on Oppenheim, Urtiew, and Stern (ref. 5).) 



variation and a time-dependent analysis of the shock- 
area interaction is carried out, the wave pattern ob- 
served within the ambiguity region is Illb. 

The minimum entropy production principle failed 
to predict the valid solution because the entropy pro- 
duction was incorrectly defined. The total entropy of 
an infinitesimal element of mass is SpA dx . If we inte- 
grate between x = — oc and x = oo at a fixed time £, 
we get the total entropy in the channel. The entropy 
production in an interval of time At is, therefore, 
given by 


/ oo 

\p( At)S( At) - p(0)S(0)} A dx (19) 

-oo 


Equation (19) can be easily integrated in closed form. 
For wave pattern IIIa 3 figure 8, we get 


A 

A t 


{P3S3 - P7S7) w r + 


S5 f 
a JO 


U5 Q-5 


pdx 


+ - [pbSba5 + PiSi{w t - U4)] 
a 


( 20 ) 


Region 7 is downstream of the reflected shock, and 
w r and wt are the speeds of the reflected and trans- 
mitted shocks, respectively. The remaining integral 


in equation (20) integrates to 



exp(— 5 5 ) 
7 



X [c k > 6 - (C - u 5 + 05)*/*] 6 - 


r* 1 

C = — + uq 


( 21 ) 


With similar results for the other two wave patterns, 
we obtain figure 14. Now, the standing shock solu- 
tion links patterns Ilia and Illb without overlapping 
pattern Illb, and the latter produces the minimum 
entropy consistent with Rudinger’s time-dependent 
computations. The figure also shows that the transi- 
tion between pattern Illb and Ilia across curve c is 
discontinuous. 


If we consider the wave patterns along curves a 
and c, figures 4 and 12, we see that the patterns 
are very similar, and we can think of curve c as the 
extension of curve a into the third quadrant. The 
same can be said of curves b and d. Curve d is not a 
boundary between two different wave patterns and, 
now that the ambiguity has been resolved, it could 
be disregarded. 

Consider the fourth quadrant, Mi < 2.068 and 
a > 1. Here we find wave pattern IVa, illustrated 
in figure 15 for Mi — 1.500 and a = 1.3. The 
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salient features are a reflected shock moving into 
the subsonic flow behind the incident shock and 
an isentropic acceleration of the flow entering the 
area contraction not sufficiently strong to generate 
supersonic flow in the small chamber. If we hold a 
fixed and increase M\ , curve e is met when reaches 
the value 1.988. At this point the flow in the small 
chamber reaches sonic conditions. Curve e, thus, is 
the locus of points separating patterns Ilia and IVa. 
If a — > oo, curve e approaches asymptotically a value 
of 1.718 for 7 = 1.4. 

Quasi-One-Dimensional 
Time-Dependent Model 

In this formulation, the discontinuous area jump 
is replaced by a steep area change defined by 

A i x ) = \ { A L + Ar) ~^( a L~ a r) tanh(cr:r) (22) 

The transition from to Ap is centered about 
x = 0 and takes place in an interval approximately 
equal to 2/<r. 


The quasi-linear form of equation (23) is 

U* 4- D(U)Uj; = Q (26) 

where 


D < u > ■ w = 


0 


1 


3(7 -3)u 2 
28 — 7 ue/p 


(3 - 7 )u 
7 e/p — 3<5u 2 


0 

28 


7uj 

(27) 

{Xrnt k ) = (x 0 + 

constant, but A t k 


We introduce a discrete grid 
n Ax, + &A t k ) where Ax is 
changes from time step to time step to satisfy the 
CFL condition (ref. 10). On this grid, we obtain an 
approximation to our dependent variable U at cell 
centers x n + ^Ax using the Roe scheme (ref. 11) to 
approximate the flux derivative in equation (23). In 
the original Roe scheme, the dependent variable U is 
interpolated to the cell faces. In our implementation, 
wc first construct the characteristic differences A W± 
from 


AW f = Ljy AU 


(28) 


where 


AU + = U n+J - Un 
AU" = U„ - U„_! 


(29) 


and lj(/ is the left eigenvector matrix of D(U) eval- 
uated with U n values. The characteristic difference 
is then limited using the minmod limiter 


A W ± 


V 


minmod jAW^jr'AW^ 

3 -K 
1 - K 


(30) 


Inside the duct defined by equation (22) we solve 
the quasi-one-dimensional Euler equations in weak 
conservation form 


Uj + F x — Q (23) 

where 


u = < 

n 

1 pu i 

> F = 

pu 

p + pv? 

Q = < 

' pu(3 y 

pu 2 P > 


[pel 


. u(pe + p). 


u(pe + p)f3 

< 4 


(24) 

and 

p = 8p(2e - u 2 ) I 



where 

TO sign(2r) / sign(y) 

mmmod[2 7 y\ = < 

l sign(^) min(|z|, |y|) sign(z) = si gn(y) 

(31) 

and K is the free constant in the kappa interpolation 
of van Leer (ref. 12), which we use to interpolate A W 
to the cell faces. In this application K — 1/3. At the 
cell faces AU is reconstructed from 

AU = Lj) 1 AW (32) 

The additional work to construct the characteristic 
differences and then the conservative variables was 
required in order to capture a strong shock. Without 
this work, the algorithm produces large oscillations 
and eventually fails. The rest of the flux evaluation 
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follows the Roe scheme as described in reference 11. 


Equation (23) is integrated in time using a three- 
stage Runge-Kutta scheme. Let 


ft(U) = At fc (Q - F x ) 

then U at time level k + 1 follows from 

u(°) = u fc 

xji 1 ) = u (0) -i- l -n 

u (2) = u( 0) + ► 

xj( 3) = u (0) + n (u (2) ) 
u fe+1 = u (3) 


(33) 


(34) 


Although the scheme allows a CFL number of 2.8, 
we have used a CFL number of 1 to avoid wiggles 
at shock waves. The overall scheme is second order 
accurate away from discontinuities. 

Results 

Comparisons between the self-similar model and 
the quasi-one-dimensional time-dependent model are 
presented in this section. The integration of the lat- 
ter is done from x — —2 to x = 2. The incident 
shock is located at x = —0.5 at t = 0. For these 
cases, a - 10 and Ax = 0.02. The first case is for 
Mi = 1.500 and a = 0.5. This case is illustrated in 
figure 5. It corresponds to a type lb pattern with a 
standing shock within the area constriction. The re- 
sults from the quasi-one-dimensional time-dependent 
solution are shown in figure 16. The reflected ex- 
pansion, standing shock, and transmitted shock arc 
clearly shown in the Mach contours. In figure 17, 
the Mach number distribution at t = 2.5 is com- 
pared with the levels predicted by the self-similar 
model. The agreement between the two models is 
good. For the second comparison, we have chosen 
conditions corresponding to figure 7, M{ = 2.500 and 
a = 0.5. At these conditions, no wave is reflected 
and a secondary shock running downstream appears. 
The expected features are clearly shown in the Mach 
contours in figure 18. The Mach number distribu- 
tion at t — 1 is compared to the self-similar solution 
in figure 19. The agreement is good except for the 
slip line in the quasi-one-dimensional time-dependent 
solution. The slip line is spread over several mesh 
points. This is a typical problem of shock capturing 
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Figure 16. Mach number contours from solution to quasi-one- 
dimensional time-dependent equations for = 1.500 and 
a = 0.5. 


Mach 

number 



Figure 17. Comparison of self-similar and quasi-one- 
dimensional time-dependent solutions for = 1.500 and 
a = 0.5. 

schemes. The third case chosen corresponds to fig- 
ure 8, M{ — 3.500 and a = 1.3. This case consists 
of a reflected shock and a rarefaction wave running 
downstream. Figure 20 shows the formation of the 
reflected shock as the left-running characteristics co- 
alesce and the formation of the rarefaction fan from 
the other family of characteristics. A Mach number 
cut at t = 1.4 is shown in figure 21. The agreement 
between the two models is good, but the compression 
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Figure 18. Mach number contours from solution to quasi-one- Figure 20. Mach number contours from solution to quasi-one- 
dimensional time-dependent equations for M< = 2.500 and dimensional time-dependent equations for M, = 3.500 and 

a = 0.5. a = 1.3. 



Figure 19. Comparison of self-similar and quasi-one- 
dimensional time-dependent solutions for M, = 2.500 and 
a — 0.5. 

behind the transmitted shock is slightly underpre- 
dicted by the quasi-one-dimensional time-dependent 
solution. For the last case, we have chosen condi- 
tions within the region of ambiguity, M* = 3.500 and 
a = 1.1. As predicted by the principle of minimum 
entropy production, the wave pattern corresponds to 
pattern Illb with a rarefaction wave running down- 



x 


Figure 21. Comparison of self-similar and quasi-one- 
dimension al time-dependent solutions for Mi = 3.500 and 


stream. The results are shown in figures 22 and 23. 
Figure 23 shows that the isentropic recompression 
produced by the area contraction is properly pre- 
dicted by the quasi-one-dimensional time-dependent 
model; however, the expansion running downstream 
shows a wiggle near its head. 
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Figure 22. Mach number contours from solution to quasi-one- 
dimensional time-dependent equations for M l = 3.500 and 
a = 1.1. 


a configuration known to be unstable. The pat- 
tern predicted in this region by numerical solutions 
of the quasi-one-dimensional time-dependent Euler 
equations is in agreement with earlier results. The 
entropy produced by the wave system was defined. 
It was then shown that the admissible pattern in 
the ambiguous region is in agreement with the pre- 
dictions of the minimum entropy production princi- 
ple. This resolved some criticisms of this principle, 
when applied to this problem, raised by Rudinger. In 
general, good quantitative agreement was observed 
between the self-similar model and the quasi-one- 
dimensional time-dependent model. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
June 18, 1991 
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